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Abstract 

A molecular model is proposed of a bilayer consisting of fully saturated DPPC and mono unsat- 
urated DOPC. The model not only encompasses the constant density within the hydrophobic core 
of the bilayer, but also the tendency of chain segments to align. It is solved within self-consistent 
field theory. A model bilayer of DPPC undergoes a main chain transition to a gel phase, while a 
bilayer of DOPC does not do so above zero degrees centigrade because of the double bond which 
disrupts order. We examine structural and thermodynamic properties of these membranes and 
find our results in reasonable accord with experiment. In particular, order-parameter profiles are 
in good agreement with NMR experiments. A phase diagram is obtained for mixtures of these 
lipids in a membrane at zero tension. The system undergoes phase separation below the main- 
chain transition temperature of the saturated lipid. Extensions to the ternary DPPC, DOPC, and 
cholesterol system are outlined. 



1 



I. INTRODUCTION 



The hypothesis that the hpids in the plasma membrane are distributed inhomogeneously 
has received enormous attention. Small domains, known as "rafts", rich in saturated lipids 
and cholesterol, have been implicated in many biological processes including endocytosis, 
transcription and transduction processes, and viral infection. The size and nature of such 
domains in biological membranes is currently unclearA. The situation in model membranes, 
however, is more transparent. Experimental studies in giant unilamellar vesicles^i^ con- 
taining mixtures of saturated and unsaturated lipids and cholesterol, which mimic the com- 
ponents of the outer leaflet of the plasma membrane, show the existence of at least three 
phases; a saturated lipid-rich gel phase, a saturated lipid and cholesterol-rich liquid phase, 
and an unsaturated lipid-rich liquid phase. The two liquid phases coexist in some regions 
of the phase diagram. The sensitivity of the phase separation to the components and their 
composition is demonstrated by the fact that mixtures of unsaturated lipids and choles- 
terol, mimicking the components of the inner leaflet of the plasma membrane, do not phase 
separate^. 

Theoretical consideration of mixtures of lipids and cholesterol have been carried out on 
binary mixtures of saturated lipids and cholesterol using models in which the enormous 
numbers of degrees of freedom of the lipids are severely reduced. In spite of this simplifica- 
tion, these models have demonstrated the importance of the preferential interaction between 
cholesterol and saturated lipids in bringing about a "liquid ordered" phase in which the tails 
are relatively well-orderedS*^. 

It is our contention that two conditions are sufficient for bringing about coexistence 
between a cholesterol and saturated lipid-rich liquid ordered phase and an unsaturated lipid- 
rich "disordered liquid" in which the tails are typically disordered. The first of these is 
that the system be below the main-chain transition of the saturated lipid to its gel phase, 
a transition driven by the ordering of the lipid tails. This transition is first-order, and is 
characterized by discontinuities in those thermodynamic quantities which are first derivatives 
of the free energy with respect to its arguments; e.g. the entropy, the area per head group, 
etc. In particular, upon introduction of an unsaturated lipid, discontinuities arise in the areal 
densities of the two lipids, first derivatives of the free energy with respect to the chemical 
potentials. In other words, phase separation between a saturated lipid-rich gel phase and 
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an unsaturated lipid-rich disordered liquid phase is a simple consequence of the main-chain 
transition. The difference in component concentrations in the two phases is expected to be 
large because the presence of the cis double bond causes the unsaturated chains to pack 
poorly with the saturated chains that are well-ordered in the gel phase. It also follows 
from the first-order nature of the main chain transition that the addition of cholesterol 
to the saturated lipid will result in a phase separation between a saturated lipid-rich gel 
phase and a liquid phase of lipid and cholesterol. Whether this liquid phase will coexist 
with the disordered liquid phase or simply transform smoothly to it as the composition of a 
ternary mixture changes depends on how different the liquids are. Thus the second sufficient 
condition for liquid-liquid phase separation is, as noted above, that the cholesterol favor 
increased ordering of the saturated lipid tails, so that the disordered tails of the unsaturated 
lipid will also tend to be expelled from the ordered fluid causing coexistence between it and 
an unsaturated-rich disordered fluid. 

In this paper, we consider a mixture of saturated and unsaturated lipids within a bilayer 
and show, as contended above, that below the main chain transition of the saturated lipid, 
there is a demixing phase separation. We do this in a model in which the lipid tails are 
treated microscopically according to the rotational isomeric states model of Flory*^ which 
permits each CH2 group to be in one of three configurations; the lowest energy trans, or 
thermally excited gauche plus or gauche minus. In the latter, the chain configuration exhibits 
a kink. Due to thermal excitation and de-excitation, this kink is transient in contrast to 
that in the unsaturated chain which is permanent. To model much of the effect of the 
interactions between the lipid segments, we follow early worb^^ and constrain the density 
of the bilayer interior to be constant, equal to a value typical of its oily components. As 
shown by these authors, imposition of this packing constraint causes the system to exhibit 
a single disordered liquid phase, characterized by a certain average number of gauche bonds 
per chain, and a certain area per chain. Larger areas could be attained by increasing the 
average number of gauche bonds in the chains, but this increases the energy and free energy 
of the system. Smaller areas could be attained by eliminating even more gauche bonds, but 
this decreases the chain entropy to the point that the free energy again increases. Thus the 
observed liquid phase occurs at an areal density at which an incremental decrease in entropy 
of the chains is just compensated by an incremental reduction in energy of gauche bonds^i. 

The constant density constraint is not adequate to bring about the main-chain transition 



3 



to a gel phase, a transition which, we have argued, is one of the sufficient conditions for 
raft formation. In order to bring about both hquid and gel phases, we not only impose 
the constant density constraint, but also include a local orientational interaction between 
normals to the CH2 planes of adjacent chains, one which tends to align them. We expect 
that inclusion of this explicit orientational interaction will capture the aligning tendency 
of the inter-chain packing. Indeed we find that this interaction does bring about a second 
minimum in the free energy, one characterized by fewer gauche bonds than in the liquid 
phase, and which we identify therefore, with the gel phase. In contrast to the liquid phase, 
an incremental reduction in entropy of the more ordered chain conformations is compensated 
here by an incremental decrease in orientational energy. 

We first consider a single component system of saturated tails consisting of fifteen 
monomeric segments, modeling those of dipalmitoylphosphatidylcholine, (DPPC). We find a 
first-order main chain transition at a temperature, T^, and calculate the degree of segment 
order as a function of position down the chain, and find good agreement with experiment. 
There is a concomitant change in area per head group as the chains lengthen in the gel 
phase. Next, we apply our technique to a bilayer of unsaturated lipids. The unsaturated 
chains we consider arc seventeen monomeric units long, with a cis- double bond located be- 
tween monomers eight and nine. This resembles the chain structure of dioleoylphosphatidyl- 
choline, (DOPC). We find no transition for bilayers composed only of this unsaturated hpid 
above the melting point of water. Finally, we consider a mixture of these saturated and 
unsaturated chains. At temperatures below T^, the system phase separates into a saturated 
lipid-rich gel phase and an unsaturated lipid-rich disordered liquid phase. Although there 
is little experimental data on this mixture, we find acceptable agreement with that which 
exists. 

This paper is organized as follows. In the next section, we introduce the theory for a 
bilayer slab composed of a mixture of two lipid tails, one saturated and one unsaturated, 
and also introduce the interaction we employ for this system. With these definitions, the 
partition function for the full mixture is specified. It is then calculated within the framework 
of self-consistent field theory in Appendix 1. Central to this theory are two self-consistent 
equations which specify two local, effective fields. The first enforces the constant-density 
constraint, and the second tends to align the local chain segment orientation. This subsection 
is written for a general system of two lipid species, one with identical saturated tails, the other 
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with identical unsaturated tails. In the following subsection, we introduce the microscopic 
model we use to describe the lipid tails. 

In section III, the general theory is applied to two specific lipids, DOPC and DPPC 
We first consider two different pure bilayers, those composed of only one or the other of 
these lipids. We present results pertinent to the main-chain transition in the DPPC bilayer 
and some geometric, conformational, and thermodynamic statistics of these membranes 
calculated in the context of our model. Of note are the deuterium order parameters for 
these lipids since they indicate the acyl chain's average conformation in the bilayer due to 
packing with neighboring lipids. In the final subsection, we present the phase diagram for the 
system of mixed lipid bilayers and some conformational statistics calculated in this model. 

II. THE MODEL AND ITS SELF-CONSISTENT FIELD SOLUTION 
A. Theory 

We consider a bilayer membrane of area A composed of A^^ saturated and unsaturated 
lipids whose tails occupy a volume V, the total volume of the hydrophobic core of both bilayer 
leaflets. Head groups are placed symmetrically, facing the external, aqueous environment 
of solvent. Temporarily, we ignore all interactions between head groups and the solvent 
molecules focusing only on the hydrophobic core of the bilayer. Acyl tails are tethered to 
the glycerol backbone of the lipid head groups on both sides of the bilayer, and the first 
monomer segments extending from the backbone define the interfaces of the hydrophobic 
core. The full width of this hydrophobic core, including both leafiets, is L. The system of 
acyl tails comprising the hydrophobic core is treated as incompressible. This assumption 
constrains the average density of monomers to a constant liquid-like value throughout the 
bilayer. The total volume is simply the sum of all monomer volumes. 

In the following, we limit our attention to lipids consisting of two identical tails. We 
further simplify the system to that of a mixture of the individual, interacting chains rather 
than a system of chain pairs extending from a common headgroup. This should make little 
difference to the chain conformations because the chains are so tightly packec^^. 

The saturated lipid has the tail structure, P—{C H2)ns-i~C H3, which we model as a chain 
of Hs linked segments extending from the glycerol backbone labeled P. All monomers except 
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the last are CH2 segments with monomeric volumes Vs{k) = = 28A^, = 1, ra^ — 1. The 
last monomer is a Cifs segment, which is assigned a volume Vs{ns) = 2i/cr^'^^^- Other than 
these volumes, the monomers are without structure, and spaced a bond length of Iq = 1.53A 
apart. The microscopic monomer number density of the saturated lipid, averaged over the 
plane of the bilayer, is 
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k=l 



(1) 

(2) 



where Zk^^{a^) is the z coordinate of the fc'th monomer of the 7'th saturated chain when 
the chain is in the configuration a^. The circumflex denotes that a quantity depends upon 
the lipid configuration. Here the single chain density, 0s^^, depends upon the single chain 
configuration a^, and the total chain density, $s depends upon the configuration S of all the 
lipids. For typographical convenience, we will simply write ^s{z) henceforth, suppressing 
the explicit dependence on the lipid configurations and allowing the circumflex to remind 
the reader of this dependence. We average over the plane of the bilayer as we shall be 
considering phases which are translationally invariant in this plane. 

The unsaturated lipid tail is similarly constructed of riu linked monomers. It has the tail 
structure P — {CH2)x ~ CH = CH — {CH2)y — CH^, which is a monounsaturated chain of 
nu = x + 2 + y + l segments. Each CH2 and CH^ segment is assigned a volume identical to 
that in the saturated lipid; each half of the —CH = CH— segment of the czs-unsaturated 
bond is assigned a volume of uqh = O.Sz^q^^. The microscopic monomer number density of 
the unsaturated lipid tails is, then, 

1 Nu 



A 



(3) 
(4) 



with volume weights. 



^0 if j < a; or x + 2 < j < riu 
0.8z/o if j = a; + l,a; + 2 
2z/o if j = nu. 



(5) 
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With these monomer volumes, the total hydrophobic volumes of single DMPC, DPPC, 
and DOPC lipids are 784A^, 896A^ and 985. 6A^, respectively, which agree well with the 
experimental values of 782A^, 913A^ and 984A^ for the fluid phases^^. We assume these 
monomeric volumes have no dependence on temperature. 

Much of the effect of the long-range van der Waals attraction and hard core repulsion 
is taken into account by the incompressibility constraint. This is clear from the agreement 
between results obtained with this approximation and molecular dynamics simulation for 
the liquid phases^J^. However an approximation based on density alone cannot capture the 
enhanced alignment of the tails which drives the main chain transition. For this we will 
need to know the local alignment of the chains. Our presumption is that the microscopic 
interactions within the system favor those configurations in which the tails are more aligned 
with one another because in such configurations the chains are more efficiently packed. We 
shall include a pairwise interaction £ between chains that reflects this assumption. 

The local orientation of the chains is conveniently specified by the normal to the local 
CH2 group. In particular, the normal to the plane determined by the /c'th CH2 group on 
chain 7 can be written 

The position vector of the anchoring carbon of chain 7 is denoted ro,^. Subsequently we 
shall refer to the ,y simply as "normals", or "orientation vectors". 

It is now convenient to introduce a local density of these normalsi^ which, for the saturated 
lipid is, 

i!s{z, u) = — ^ ^ ys{k)6{z - Zk^^)5{u - Ufe,^). (7) 

^ 7 k=l 

Again the circumflex is a reminder that this density depends explicitly on the chain co- 
ordinates Zfc,7 and normals ,y in a particular configuration of lipids. There is a similar 
expression, for the local density of normals coming from unsaturated chains. 

A local pairwise interaction between tail segments which depends upon their orientations 
can be written, 

^JdzJ dndW[^,{z,u) + ^u{z,u)]V{u,u')[^s{z,u) + 4!u{z,n')]. (8) 
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Here, / du denotes an integration over the full solid angle. It is important to note that all 
interactions are local, i.e. between bonds, and are of the same form and strength irrespective 
of the length or degree of unsaturation of the chains to which the bonds belong. 

Because the normal to the plane of the bilayer, c, singles out a spatial direction, we 
expect the interaction to depend upon the angles that the segments make with this normal 
V{u,u') — \^(u',u) V{u • c, u' • c). Irrespective of the particular form chosen for this 
interaction, within the mean-field approximation which we introduce shortly, the dependence 
of the free energy of a single chain on the normal can only arise through functions of u • c. It 
follows from this and the fact that the interaction energy, E, is pairwise that the free energy 
of the system of chains can only depend upon a product of such functions. Anticipating 
this, we take for the interaction a simple form which already exhibits this product nature, 

V(u • c,u' • c) -Jf(n ■ c)/(u' • c) = - Jf{cose)f(cose'), (9) 

where f{cos9) is normalized. Further, in the absence of headgroup interactions which we 
have ignored, the chains will align normal to the plane of the bilayer which implies that 
/(cos 9) be chosen a monotonically decreasing function of its argument. With this orienta- 
tional interaction, the lowest energy configuration of two neighboring orientational vectors 
is one in which they avoid the penalty of hard-core interaction by both aligning with the 
bilayer normal. 

With this choice of orientational interaction, it is convenient to define 

E,{z) = J d^sin^^',(2;,u)/(u-c) 

= tEU^), (10) 

^7=1 

ris-l 

= ^s{k)6{z - Zk,^)f{uk,j ■ c), (11) 

k=l 

with a similar expression for the unsaturated lipid tails, so that the interaction energy can 
be written in the transparent form, 

^--i^Jd^[^s{z) + E^{z)f. (12) 

For /(cos 6') we choose the simple functional form 

/(cos^) = ^^^(cos^e)"^ 

mexp(— m^^), (13) 
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which is approximately a Gaussian with width m~2. The prefactor is a normahzation 
coefficient. If the interaction, Eq. (|T2|l . were not local but were an average over the lipid 
chain, and if m were unity in the expression above, then this interaction would be comparable 
to the leading order term of the interaction Marceljaii^ considered for saturated lipids. The 
local form, and higher power dependence of the alignment with the bilayer normal was first 
considered by GrueniS. 

Two parameters determine this interaction between segments, its strength J and its 
angular range. The former is set by requiring that the main-chain transition temperature 
for the chains under consideration agree with experiment. The latter may be estimated by 
considering the average angular alignment of an acyl tail in the bilayer. Crudely, a lipid in 
the bilayer is, on average, confined to a cylindrical shape. A simple estimate of its alignment 
with the normal is then 6*0 ~ where d is an average cross-sectional area for an acyl 

chain and / is an average length in a typical membrane leaf. Comparison of the Gaussian 
width to this estimate, leads to m ~ 14 — 18. We use m = 18 in our calculations. 

Now that the interactions between chains are defined, the Helmholtz free energy of the 
system, subject to the incompressibility constraint of constant density, can be obtained 
within self-consistent, or mean field, theory. This is done in Appendix I with the result 
l3Fr,ft{T,N,,N^,A) _ , ^ ^ n , f^-J 



^ ps In - P« In + — j dz[ps < ^siz) >, +p„ < ^^iz) >„] 
H In -— H m 
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2 psEkJ^sik) + puEkJ^uik) 2 psllk^s{k) + puT.k^u{k) 
^ ^ dz 7r(z), (14) 



where ps = N^/A and pu = N^/A are the areal densities of saturated and unsaturated chains. 
Here < ^s{z) >s is the average in the ensemble of a single saturated chain 

<Uz) >s=Tr{^yPsUz) (15) 

where the trace is over all configurations of a single saturated chain, and the probability 
distribution Pg is given by 

= -^exp(-/3^i-- f[4)s{z)7r{z)+Uz)b{z)]dz\ (16) 

with Qs the single saturated chain partition function 

Q,(7r,6) = Tr|,|exp|-/3ifi - \^s{z)t^{z) + Uz)h{z)]dz^ , (17) 



and Hi the intra-chain part of the Hamiltonian. The self-consistent equations which deter- 
mine the fields ni^z) and h{z) are 

I = Ps < (t)s{z) >s +Pu < (i)u{z) >n, (18) 

b{z) = -(3J[ps < L{z) >s +Pu < L{z) >u]- (19) 

The first of these enforces the incompressibility constraint. 

Note that, because of the incompressibility constraint, the volume is not an independent 
variable so that the Helmholtz free energy, Eq. ()14|) . depends only on the temperature, the 
numbers of each lipid, and the area. The interpretation of the free energy is straightforward. 
The first two terms are simply the single chain Gibbs free energies multiplied by the number 
of chains per unit area. As is well known, these terms double-count the pairwise interaction 
energy, and this is corrected by the third term. The next two terms are the entropy of 
mixing, and the last is simply —(3pV/A which relates the Gibbs free energy of the preceding 
terms to a Helmholtz free energy. 

Other than excluding tail configurations protruding into the solvent region, the effects 
of the aqueous medium and head groups have thus far been neglected. It is evident that 
the interfacial effects can also be somewhat complex. For example, tilted solid phases are 
often attributed to a mismatch between a smaller average areal footprint of the packed tails 
and a larger headgroup. On the other hand, a small headgroup mismatched with larger, 
fluid-like tails attached to it can yield costly interfacial interactions between the solvent and 
the hydrophobic core, resulting in a large interfacial tension. These effects are beyond the 
scope of our calculation. However, so as not to overlook the contributions from the head 
groups completely, we choose a minimal model for them. 

Ignoring the molecular details of the head groups and their interactions with the aqueous 
environment, we focus on two physical attributes that result from these interactions. The 
first is an interfacial energy per unit area, which we include as a simple phenomenological 
form, JqA. Thus the total free energy per unit area is 

/^(T,p„pJ = F^^t/A + 7o, (20) 

where F^ft is given in Eq. ()14|) . The value of 70 is taken to be on the order of the energy 
per unit area of an oil-water interface. We use 70 = 0.12A;bT/A^ in our calculations, which 
is roughly the experimental value. The second effect attributable to the head groups is the 
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rigidity they provide the acyl taiL The hpid tails are tethered to the glycerol backbone, 
providing a fixed endpoint for the chain. The head groups tend to order the first segments 
of the lipid tethered there, an effect which we model by a surface field that couples to the 
first bond of the lipid. To model this, we include in Hi a surface field term r7cos(ui ■ c). 
The magnitude of the surface field rj is tuned so that the average angle of the first bond is 
approximately equal to that obtained from NMR measurements for the lipid. 

B. Microscopic Model 

To evaluate the single lipid partition function in external fields, Eq. (fTTj) above, we specify 
a microscopic model suitable for describing the noninteracting single chain configurations 
and their energy. Hi. To this end, we employ Flory's Rotational Isomeric States (RIS) 
model^ for the lipid acyl chains. This model places a subsequent monomer in the acyl chain 
at one of three discrete possible orientations from its ancestor segment; the trans, gauche 
plus, and gauche minus orientations. According to this model, the gauche bonds are local 
excitations that contort the chain and increase the energy by tg = 500 cal/mol. Thus Hi 
of a given configuration is simply equal to eg multiplied by the number of gauche bonds 
in that configuration, and also includes the contribution from the surface field described 
above. The lipid conformation with all trans bonds is a highly stretched chain, whereas a 
lipid tail with several gauche bonds is coiled. Given the location of the first bond in the 
acyl chain, all positions of the subsequent — 1 generations are determined by this model, 
which locates segments in a tetrahedral fashion, with skeletal angles of 112°, and dihedral 
angles of {0°, ±120°} for trans and gauchei orientations respectively. Without restrictions 
on this set, there are 3""^^ configurations for each fixed, first bond orientation. 

The CIS- double bond in the unsaturated chains slightly modifies this model. This bond 
is a quenched defect in the chain, effectively placing a kink in the tail about the location 
of the double bond. The energies and configurations of the segments adjoining the double 
bond as well as the —C = C— bond length differ slightly from the saturated links. For 
the segments associated with the double bond, we employ the RIS configurational model of 
c2s-l,4-polybutadiene^^, which possesses the same atomic structure as the double bonds in 
these lipid chains. 

We also impose several restrictions on the lipid configurations that reduce the number 
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of allowed ones significantly. First of all, configurations that protrude into the aqueous en- 
vironment are discarded. A configuration with any two, non- consecutive monomer centers 
less than a bond width apart is considered self-intersecting and also discarded. Configura- 
tions with consecutive gauche minus (plus), gauche plus (minus) bond sequences have large 
hydrogen- hydrogen steric overlaps^, a phenomena referred to as the "pentane effect." These 
configurations are also neglected. Each single lipid partition function is then evaluated by 
enumerating all allowed configurations in the RIS model, from many different first bond ori- 
entations and locations. It should be noted that two different first segment orientations will 
have different numbers of configurations stemming from them since they will have different 
number of configurations that penetrate the aqueous environment. 

We specify in Appendix II how the probability distributions Pg and P„ and the single 
lipid partition functions, Eq. (fT7j) . are calculated within the context of the RIS model and 
the configurational restrictions. We also discuss there a few other details of the calculation 
of the free energy. 

III. RESULTS 

A. Single component DPPC and DOPC bilayers 

We first apply our method to a tensionless bilayer composed only of DPPC lipids. We 
solve the self-consistent equations for a given temperature and under the condition that 
the surface tension be zero. The high temperature liquid phase undergoes a main chain 
transition to a low temperature gel phase, and by adjusting the interaction strength J, we 
set this transition temperature at Tm ~ 315 K, essentially the experimental one^^. The value 
of J remains fixed at this value throughout all subsequent calculations. In Fig. 1 we show 
the Helmholtz free energy per two-chain lipid 



in units of fc^T plotted vs. the area per two-chain lipid, 2A/Ns = 2as, at the temperature 
of the main chain transition. It is straight forward to show that 




(21) 
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{S/Ns)dT], 



(22) 
(23) 
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with S the total entropy. The two minima in the figure occur at the areal densities of the 
two coexisting phases. Because these states are minima of the free energy per hpid, their 
surface tensions vanish, as seen from Eq. ()23p. Further as they have the same free energy 
per hpid at zero surface tension, they have the same chemical potential, as seen from Eq. 
fl22|) . Thus these phases are indeed at coexistence. The main chain transition is clearly seen 
to be first order, marked by a discontinuous jump in the area per lipid. Table H] summarizes 
a few geometric aspects of each state and of the transition. We find satisfactory agreement 
with experiment for the geometric properties. The fiuid phase is characterized by a greater 
area per lipid than is the gel phase. The hydrophobic width of the bilayer is thinner, filled 
with coiled lipid tails with, on average, somewhat more than four gauche bonds per tail. 
The gel state consists of more ordered lipid chains that are tightly packed, having a smaller 
area per lipid and fewer gauche bonds. 

One thermodynamic quantity that is sensitive to the packing in the hydrophobic core 
is the deuterium "S'cd" order parameter profile, which is commonly measured in the fiuid 
state by nuclear magnetic resonance (NMR) experiments. The calculated deuterium order 
parameter for our model DPPC bilayer is shown in Fig. 2 for both fiuid and gel states. The 
fiuid phase is compared to experimental data from reference 20. Collective organization of 
the lipid tails during the transition from the fiuid to gel state is indicated by the large jump 
in chain alignment. Experiment and simulation distinguish sn — 1 and sn — 2 lipid tails, 
which we do not, so the calculated profiles in Fig. 2 should be considered as an average over 
both acyl tails. Comparing to experiment, we see that the calculated fiuid phase is slightly 
more disordered than the experimental one. 

Applying our calculational method to a bilayer consisting only of mono-unsaturated, 
DOPC model lipids, we find a fiuid phase and no transition above the freezing point of 
water. This result is consistent with experiment, which indicates a transition temperature 
for DOPC of Tm ~ 250 K. The calculated fiuid state has an area per lipid of 70. 4A^ and a 
hydrophobic width of 28. OA. This lipid area is somewhat smaller than the experimental value 
of 72. 5A^ with also a thinner hydrophobic width of 27. lA. The calculated deuterium order 
parameter is shown in Fig. 2. It should be noted that these profiles are calculated for a lower 
temperature, T = 305 K, than the saturated lipids. The order parameter displays a notch 
in the profile about the cis- double bond, one which is also seen in MD simulation^SiSLS^. 
Also, there is a slight even-odd effect in the order parameter which we do not observe in the 
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saturated chains. 



B. Composite bilayers of DPPC and DOPC 

From the calculation of the free energy, Eq. ()14|) . the phase behavior for a mixture of 
both lipids, DPPC and DOPC, can be determined. Phase coexistence is determined by 
examining the auxiliary function, 

7 = min {/a- Kps - Kpu} ■ (24) 

ps ■jpu 

When the and A„ are adjusted so that 7 exhibits two minima with the same value of 
7, then there is two-phase coexistence. At this point, the values of A^, \u are equal to the 
values of the chemical potentials of the saturated and of the unsaturated lipids, /i^ and /i^j, 
in the coexisting phases and 7 is equal to the surface tension, 7, of each phase. We further 
require that this surface tension be zero in each phase as there is no constraint on the bilayer 
area. The four requirements that the chemical potentials and surface tension be equal in 
each phase and that the tension be zero determines the two areal densities, ps and pu-, in 
each of the two phases. Note that the bare surface tension, 70, causes the system to contract 
to the point at which the tension due to packing just counteracts the bare tension leaving a 
net tension of zero. 

The calculated phase diagram of the system is shown in Fig. 3. For small concentrations 
of DOPC, we find a gel state rich in DPPC below the main-chain transition temperature. 
A region of fluid/gel coexistence is found which is rather wide because of the poor packing 
of the unsaturated lipids. These results are in reasonable agreement with experiment^. 

The deuterium order parameters for the two coexisting phases at T = 293 K are shown 
in Fig. 4. The gel phase shows significant ordering of both lipids, given by the large 
absolute magnitude of both order parameters. The disordering effect of the double bond 
is more prominent in the ordered phase, with affected bonds showing a large notch. This 
is noteworthy because it implies that the tight packing about the unsaturated lipids in 
the ordered phase tends to enhance the disorder about the unsaturated bond. The order 
parameters in the liquid- crystalline phases are also shown in Fig. 4. 
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IV. SUMMARY 



We have formulated a microscopic model of a bilayer consisting of lipids, one which 
incorporates two important effects of the interactions between chains; the tendency of the 
chains to create an incompressible hydrophobic interior, and their cooperative tendency to 
order. The solution of the model obtained from self- consistent field theory reproduces many 
salient features of the lipid bilayer. The system of chains which mimic DPPC undergoes a 
first-order main chain transition from a disordered fluid state to a more ordered gel phase, 
a transition characterized by a reduction in the number of gauche bonds and the area per 
chain. The ordering of the chains in the liquid phase as determined by the deuterium 
order parameter is in rather good agreement with experiment. We have also calculated the 
deuterium order parameter in the gel phase. 

We then applied the model to a system mimicking the mono-unsaturated chains of DOPC. 
Utilizing the same interaction as before, we find that this system does not undergo a main 
chain transition above zero degrees centigrade, a result in agreement with experiment. The 
lowering of the main-chain transition is due to the double bond in the chain which disrupts 
the tendency to order. The deuterium order parameter profiles for this system show the 
characteristic reduction at the position of the double bond. 

Lastly wc examined a mixture of saturated and unsaturated lipids and showed that the 
system phase separates into a saturated-lipid-rich gel phase and an unsaturated-lipid-rich 
liquid phase at temperatures below the main-chain transition temperature of the saturated 
lipid. The separation arises from the first-order main chain transition. The large difference 
in areal densities is due to the lipid architecture; the unsaturated lipids are rejected from 
the gel phase because they do not pack into it very well. 

Several pieces of physics have been ignored in our calculation, such as the explicit inter- 
action between head groups and, of course, the effect of fluctuations. We expect that such 
effects will alter the specific values of thermodynamic quantities, but not the main phenom- 
ena. In particular the phase separation into a saturated-rich gel and an unsaturated -rich 
liquid phase will remain below the main chain transition of the saturated lipid. As noted 
earlier, this is important as this behavior is one of two ingredients for a sufficient theory of 
liquid-liquid coexistence in a ternary mixture which includes cholesterol. We are currently 
investigating the role of this third component within the framework of our model. 
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VI. APPENDIX I. SELF-CONSISTENT FIELD THEORY 



The exact Helmholtz free energy of the system is 



H + In Pex\ (25) 

where the trace is a sum over all distinguishable configurations of all chains, two chains per 
molecule, and Pex is the exact probability distribution function for the system. It satisfies 

Tr{j:yPe, = 1. (26) 

The Hamiltonian of the system is 

H = Y: Hi,, + E Hi,, - — / dz[- E i,,{z) + ^ E L,,{z)r, (27) 

7=1 r]=l J ^^=1 A 

where H\^.^ contains the intra-chain contributions to the energy of chain 7, such as those 
arising from gauche bonds and the coupling to the surface field. Self-consistent, or mean field, 
theory consists of approximating Pex by the product of one chain probability distributions. 



-iy\-Yy-ViViPs,,Pu,, (28) 

^ ^ ^ ^ / 7=lr?=l 

where the probability distribution of the single saturated chain is written 

n,7 = TT I -/^^I'T - E ^ [n(^fe,7) + 5.(^fe,7)/(ufe,7 • c)] I , (29) 

= ^exp|-/5^i - \^,{z)^{z) +Uz)B,{z)\d^ . (30) 

with the fields 11 and Bg to be determined, and Qs the single saturated chain partition 
function 

a(n, B,) = rr{,} exp - E ^ [n(^fe) + Bs(zk)f(uk • c)]| , (31) 

V k ^0 J 

= Tr^^y exp - [M^Mz) + Uz)B,{z)]dz'^ . (32) 
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Here 

Us 

4>s{z) = '^Us{k)6{z - Zk), (33) 
fc=i 

defined previously in Eq. the trace is over all configurations of a single saturated chain, 
and the index specifying the particular chain has been dropped. The expressions for P„ and 
Qu are essentially the same except that the segment volumes are those of the unsaturated 
chain, and the trace is over its configurations. Note that because 

Tr{,}P, = 1, (34) 

and similarly for and because the sum in the normalization of the exact probability 
distribution, Eq. (j26|) . is over all distinct states, the factor of (A^s/2)!(iV„/2)! is needed in 
Eq. ()28|) so that the normalization will be satisfied within the mean field approximation. 

The Helmholtz free energy of Eq. ()25|) can now be obtained directly within the mean 
field approximation of Eq. ()28|) . One obtains 



f PJ ^ ^ ^ 



P£ _ _]_ f,(3J 
A ~ uo 

+ [ps < (ps >s +Pu < <Pu >u]^{z)]dz 



- [ps\n Qs + Pu\n Q„] + yln iV, + ^ln iV„ (35) 
where ps = Ns/A and p„ = Nu/A are the areal densities of saturated and unsaturated chains, 

< (f)s >s= Tr{a}Ps(j)s, (36) 

and similarly for the other single-chain ensemble averages. Again chain indices have been 
dropped. 

The unknown fields 11(2;), Bs{z), and Bu{z) are obtained from three conditions. The first 
is that the interior of the bilayer be incompressible which implies 

ps < 4>s{z) >s +Pu < 4>u{z) >u= 1- (37) 

The second and third are that the free energy be an extremum with respect to the undeter- 
mined functions < C,s >s and < >„. This leads to the equations 

Bsiz) = B^{z) (38) 
= (3J[ps < L{z) >s +Pu < L{.z) >u]- (39) 
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We shall denote the functions which satisfy the self-consistent Eqs. ^7\i to (jSHl) as 7i{z) and 
b{z) respectively. 

Substituting the fields 7r{z) and b{z) into the free energy, one obtains 



l3F 13 J f " 

— = -pslnQs -p„ln(5„ + — J dz[ps < is{z) >s +Pu < Liz) >u 



2 



2uo 

+ ^lniV, + ^lniV„-l / dz7c{z). (40) 

The first two terms are simply the single chain Gibbs free energies multiplied by the number 
of chains per unit area. (They are Gibbs free energies as the pressure-like vr is the independent 
variable.) These terms double-count the pairwise interaction energy, and this is corrected 
by the third term. The fourth and fifth terms are the entropy of mixing. Thus the first 
five terms are the mean field approximation for the Gibbs free energy per unit area of the 
system. The last term is simply —[3pV/A, with p the pressure, and relates the Gibbs free 
energy to the Helmholtz free energy. 

The contribution of the entropy of mixing can be put in a more convenient form as follows. 
We note 

(AT, Y: m + N^Y1 ^(k))/V = 1, (41) 



so that 

where ln(V/z/o) is constant within the canonical ensemble. We multiply this constant by 
{Ns + Nu)/2 and subtract it from the free energy above because this term, linear in Ns and 
Nu, only shifts the definitions of the chemical potentials. We therefore arrive at the final 
expression for the mean field free energy 

PF I f3J f " 

= -PslnQs -p„lnQ„ + — j dz[ps < is{z) >s +Pu < Liz) >«]^ 

^ Pll^ Pf^ + ^ In ^"'^^ 



2 PsY.k^s{,k) + puY.k^u{,k) 2 psY.k^s{,k) + puY.k^u{,k) 
^ ^ dz 7r(z). (43) 



VII. APPENDIX II. THE RIS DISTRIBUTION 



This section clarifies how the sums over single lipid configurations, which are needed 
to evaluate single lipid ensemble averages, Eq. (fT3|l. are carried out within the rotational 
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isomeric states (RIS) model, and how the calculation of thermodynamic quantities is accom- 
plished. In the RIS model®, the bond between the carbon at position r^-i and , > 2, 
can be in one of three configurations, the lowest energy trans or the gauche plus or gauche 
minus with an energy higher than that of the trans by an amount eg. There can be no gauche 
bond between the anchoring carbon and the first carbon in the chain by definition. Further 
an apparent gauche bond between the first and second carbons can actually be obtained by 
a solid body rotation of a trans configuration. Within the RIS model, the sequence of bond 
orientations down the chain, the location of the anchoring carbon, zq, and the orientation 
of the first normal, ui, completely define the configuration of the chain. 

The probability distribution, Eq. (pUj) . of a saturated chain in configuration a can now 
be written 

1 ris 

71s 

Qs = Tr{a} n ^s,fc(a) 

k=l 

Ps,k{a) = exp(-/?i^e//,fc), (44) 
where the effective Hamiltonian is defined as 

(3Heff^i = r7cos(ui ■ c) + ^^^[7r(zi) + h{zi)f{ui ■ c)], 
PHeff,2 = ^[7r(^2) + b{z2)f{u2 ■ c)], 

pHeff^k = 9k(^g + — [-^(zk) + b{zk)f{uk ■ c)], 2 < < - 1 

PHeff,ns = gus^g + -T^i^nJ, (45) 

^0 

where t] is the surface field introduced earlier which couples to the first normal, and Qk 
equals unity if the bond between the carbon at Tk-i and rk is gauche, and zero otherwise. 
For unsaturated lipids, the cis-unsaturated bond alters the bond weights. These energetic 
weights and configurational properties are specified in Marki^. 

Not all possible bond configurations in a chain are acceptable. Those in which a gauche 
minus immediately follows a gauche plus or vice versa have large steric overlaps^, and are 
discarded. Further those configurations in which the chain intersects itself are also discarded 
as are those in which the chain pierces the aqueous plane and enters the region of solvent. 
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To carry out the sum over single lipid configurations, we first enumerate all configurations 
of the chains within the RIS model for a given position of the first monomer and orientation 
of the first bond. Under this constraint, the RIS configurations of these lipids are generated 
once and stored in a linked list with a tree structure. This structure allows us to store a large 
number of conformations and perform configurational sums using roughly a quarter of the 
number of operations required by a linear storage array. We then sample orientations and 
translations of this first bond by rotating it about the bilayer normal and then translating the 
chain along the bilayer normal. In particular we consider rotations of the lipid configurations 
as a solid body, sampling two Euler angles: that to the aqueous plane from the normal, 9, 
and 0, and about the body axis. 

Translations of the whole chain normal to the bilayer are chosen in a small interval 
{—Az, Az}, with Az ~ 1.5A, about the aqueous plane. We choose the first segments of the 
chains to originate at one of four evenly spaced locations in this window about the aqueous 
plane, matching a discretization of the bilayer we employ, described below. The allowed 
Euler rotations about the first bond varies with these origins. For example, the maximally 
inserted chains in this translational grid have somewhat more than the hemisphere of allowed 
rotations without penetrating the aqueous environment, whereas the maximally extracted 
ones have less than the hemisphere. 

For each of these origins in the translational window, we sample the allowed solid angles. 
Using a shorter saturated lipid, 12 or 14 segments long, as a benchmark, we may sample 
both Euler angles exhaustively. A large uniform sampling, (~ 3500 total angles), indicates 
the dense phase has a preference for two orientations, 6' ~ 0, corresponding to the first 
bond in the trans state, and to a lesser extent, 6 ~ 60° corresponding to a first bond 
gauche isomerization. Since such an exhaustive sampling of angles for longer lipids is not 
computationally feasible, we must substantially reduce the number of angles sampled. 

One way of accomplishing this is to utilize a set of angles determined from equidistant 
points on a sphere. This accomplishes a uniform density distribution on the sphere, each an- 
gle subtending nearly the same portion of the total solid angle. These angles are determined 
from a simulation of repulsive points confined to a sphere. We then increase the density 
of configurations for small 6, until we nearly reproduce the statistics from the exhaustive 
sampling above. In order to correct for the bias imposed by the selection of angles, the 
evaluation of the trace such as in Eq. (|T5|) includes a factor proportional to the amount 
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of the total allowable solid angle they subtend. The dense set of angles near ^ ~ each 
subtend a smaller amount of area and have a smaller weight. 

With these rotations and translations, and discarding self-intersecting and other restricted 
configurations, the complete computational set includes on the order of 10^ configurations 
for roughly 150 total initial bond angles for each lipid species. 

To determine ensemble averages of position dependent quantities, such as 

(psai^, Zka) = ^s{k)5{z - Zk^-y), (46) 

k=l 

of Eq. © or 

= J2 ^sik)6{z - Zk,^)f{uk,^f ■ c), (47) 

k=l 

of Eq. ()11|) . we discretize the space of the continuous variable z into 2A/'z levels zi of width 
6z = L/2N'z- A saturated chain which is without gauche defects and which is aligned 
with the bilayer normal sets a limiting, maximal half- width, L/2, for the hydrophobic core 
of a one-component lipid bilayer. Typically, we divide the hydrophobic width into slices 
5z ~ 0.75A to achieve a suitable resolution for our calculations. Free energy calculations for 
various bilayer widths are carried out by changing the number of slices, 2A/'^, but maintaining 
the slice width 5z. By changing the bilayer thickness incompressibly, we simply vary the 
average areal density of the lipids. The equilibrium states of the system are typically not 
among the discrete points we calculate, which may be seen in Fig. (1). Structural properties, 
such as the Scd order parameters, are sensitive to the areal density. To estimate the order 
parameters, shown in Figs. (2) and (4), we first estimate the areal density of the phase 
from the free energy calculation, and then vary 5z so that the bilayer acquires nearly the 
equilibrium density of the phase. In practice, changing 5z slightly affects the packing of the 
chains in the hydrophobic core so several values of 5z are used, and the order parameters 
are averaged over these various discretizations. This practice introduces an uncertainty in 
these order parameters which we estimate as ±0.02. 

The calculation is now straightforward. For given surface compositions ps and pu and 
temperature T, the densities < (psizi) >s, < 4>u{zi) >u, < is{zi) >s, and < is{zi) >s are 
evaluated in the ensemble with the probability distribution Pg of Eq. (fTB|) or Pu containing 
given fields 7i{zm) and b{zm)- These fields must satisfy the self-consistent equations, Eqs. 
fll8|) and (fT^ . which contain the densities. This set of coupled, nonlinear, equations can be 
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solved using a standard algorithm. Due to the mirror symmetry of the bilayer, it is sufficient 
to solve the self-consistent equations in the half-bilayer. To do this, chain interdigitation 
must be correctly accounted for in lipid configurations that extend past the midplane into 
the other leaf. If the configuration is long enough to interdigitate the other leaf, its monomer 
densities are refiected and the 'image' of the configuration counted in the statistics. The 
free energy, Eq. (fT^ . is then calculated from which all thermodynamic information may be 
derived. 

The parameters we have chosen for the calculations below are as follows. The saturated 
lipid chains we consider are Us = 15 monomers long, mimicking DPPC. An interaction of 
strength J/T^ = 0.09 induces the main-chain transition at its observed value of = 315 
K. The unsaturated lipids have = 17 monomers with a cis-9 double bond, between the 
eighth and ninth segments. This is the chain structure of DOPC. For both species of lipid, 
we find a surface field of rj/ksTm ~ 1-03 sufficient to obtain reasonable values of the order 
parameter for bonds nearest the aqueous plane. 
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TABLE I: Mean-Field and Experimental Results for the DPPC Bilayer. Experimental results are 
from Nagle and Tristam-Nagle^. 



Quantity 

Fluid Area/lipid T = 323 K (A^) 

Gel Area/lipid T = 293 K (A^) 

Fluid Hydrophobic width T = 323 K (A) 

Gel Hydrophobic width T = 293 K (A) 

Gauche Bonds (gel, T = 315 K) 

Gauche Bonds (fluid, T = 315 K) 



Exp't Mean Field 

64.0 67.0 
47.9 49.9 
28.5 26.7 
34.4 35.9 

2.1 
4.3 
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VIII. FIGURE CAPTIONS 



Fig. 1 Free energy per two-chain saturated lipid, in units of /c_b^, as a function of 
area per two-chain hpid. The temperature is essentially that of the main chain transition, 
approximately 315 K, at which the two phases are in coexistence. In each phase the surface 
tension, 7, vanishes. The solid line is to guide the eye. 

Fig. 2 Various calculated deuterium order parameters for a DPPC and a DOPC bi- 
layer are compared with experiment and simulation. The top figure shows the calculated 
order parameters in coexisting liquid and gel phases for a DPPC bilayer at the transition 
temperature. The middle figure compares the mean-field deuterium order parameter for the 
fluid phase at (T = 323 K) to that obtained from experimenlj^ (T = 323 K) and from MD 
simulatioii^^ (T = 325 K, sn-2 chains) for the fluid phase. The bottom figure compares our 
results for the fluid DOPC bilayer to simulation22i at = 305 K, sn-2 chains. The indices 
have been shifted from the text in accord with most experiments and simulations. 

Fig. 3 Calculated phase diagram for a bilayer composed of DPPC and DOPC as a function 
of the volume fraction of DOPC. The error bars are an estimate of the uncertainty in the 
calculation introduced by the bilayer discretization. 

Fig. 4 Deuterium order parameters for a composite bilayer of DPPC and DOPC at 
Tm = 300 K in each of the coexisting phases. The order parameters for DOPC are identifiable 
by the notch half way down the chain. The upper part of the figure shows the gel phase 
order parameters, and the lower part shows the liquid phase order parameters. The indices 
have been shifted from the text in accord with most MD simulations 
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R. Elliott et. al, Fig. 2 
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R. Elliott et. al, Fig. 3 
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